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A generalized trial wave function termed as the “multi-Dj ansatz” has been developed to study 
the ground state of the spin-boson model with simultaneous diagonal and off-diagonal coupling 
in the sub-Ohmic regime. Ground state properties including the energy and the spin polariza¬ 
tion are investigated, and the results are consistent with those from the exact diagonalization and 
density matrix renormalization group approaches for the cases involving two oscillators and two 
baths described by a continuous spectral density function. Breakdown of the rotational and parity 
symmetries along the continuous quantum phase transition separating the localized phase from the 
critical phase has been uncovered. Moreover, the phase boundary is determined accurately with the 
corresponding symmetry parameters of the rotational and parity symmetries. A critical value of the 
spectral exponent s* = 0.49(1) is predicted in the weak coupling limit, which is in agreement with 
the mean-field prediction of 1/2, but much smaller than the earlier literature estimate of 0.75(1). 
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I. INTRODUCTION 

The paradigm of a quantum spin interacting with its 
dissipative environment has drawn sustained research in¬ 
terests in a variety of fields including quantum compu¬ 
tation (Tj- 3}, spin dynamics 0-Q, quantum phase tran¬ 
sitions I7H10I] . charge transfer in biological molecules 

g d and impurity effects in magnetic materials [13- 
Among the most popular models employed in this 
regards is the spin-boson model [lgj] that describes a two- 
levcl system, i.e., a spin 1/2, coupled linearly to an en¬ 
vironment represented by a set of harmonic oscillators. 
The coupling between the system and the environment 
can be characterized by a spectral function J(w). This 
model is known to exhibit rich ground state properties. 
In particular, if the bath is characterized by a gapless 
spectral density J(w) ~ 2 aco s , a quantum phase transi¬ 
tion is expected to appear, separating a non-degenerate 
“delocalized” phase from a doubly degenerate “local¬ 
ized” phase due to the competition between tunneling 
and environment-induced dissipation. Depending on the 
value of s, there exist three distinct cases known as sub- 
Ohmic (s < 1), Ohmic (s = 1) and super-Ohmic (s > 1) 
regimes. Recent theoretical studies UE3 show that the 
transition is of second order in the sub-Ohmic regime and 
Kosterlitz-Thouless type in the Ohmic regime [/|. In the 
super-Ohmic regime, however, there is no phase transi¬ 
tion. 

A number of studies have investigated extensions of 
the standard spin-boson model, for example, to a two- 
spin system involving a common bath [17, |l8| or two 
independent baths [19j, and to a single spin coupled 
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to a bath with simultaneous diagonal and off-diagonal 
coupling [20]. We have recent studied a two-bath spin- 
boson model Hil S3, shown schematically in Fig. Oja)], 
where X and Z denote the diagonal and off-diagonal cou¬ 
pling, respectively, and the arrow represents a spin. The 
bath spectral densities can be described by J 2 (w) = 
2 aujl~ s co s , J x (lo) = 2/3w/ _s w s , where a and /3 are the 
dimensionless coupling strengths, and s and s denote the 
spectral exponents characterizing the two baths coupled 
to the spin diagonally and off-diagonally, respectively. 
Possible realizations of such two-bath model include im¬ 
purities in a magnet coupled to two spin-wave modes or 
two sources of dissipation flU [2(1 , excitonic energy trans¬ 
fer processes in natural and artificial light-harvesting sys¬ 
tems [27} , electromagnetic fluctuations of two linear cir¬ 
cuits attached to a superconducting qubit [28j - [30j ]. two 
cavity fields coupled to a SQUID-based charge qubit [3l} . 
and the process of thermal transport between two reser¬ 
voirs coupled with a molecular junction [32}. 

In the two-bath model with s = s, a = (3, studies 
based on the perturbative renormalization group theory 
predict the presence of two phases, namely, the “criti¬ 
cal phase” and the “free phase,” in the absence of bias 
and tunneling [24} (33} [H} . Very recently, the existence 
of a “localized phase” in the two-bath model was dis¬ 
covered numerically in the strong coupling regime [35} . 
The schematic of the phase diagram that emerges is 
shown in Fig. [ljb). A continuous quantum phase tran¬ 
sition separating the localized phase from the critical 
phase was claimed to exist only for the spectral expo¬ 
nent s* < s < 1, and a critical value of the spectral ex¬ 
ponent, s* = 0.75(1), was estimated from the density ma¬ 
trix renormalization group (DMRG) calculations. When 
s > 1, the impurity behaves as a free spin in the so-called 
free phase [35}. The phase boundary was determined 
from the response to the external field (i.e., the bias or 
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Figure 1: (Color online) (a) Schematic of the two-bath spin- 
boson model. A single spin is immersed in two indepen¬ 
dent baths with simultaneous diagonal coupling (Z) and off- 
diagonal coupling (X). (b) Schematic plot of the phase di¬ 
agram of the spin-boson model with two identical bosonic 
baths s = s and a = /3. Where s ( s ) and a (/3) represent the 
spectral exponents and coupling strengths, respectively, for 
the spectral density functions J z {u>) Three different 

phases (localized, critical and free) are displayed in the as 
plane with two critical values of the spectral exponents, s * 
and 1.0. 


tunneling) perpendicular to the bath plane. However, the 
localized-to-delocalized phase transition will occur under 
the external field, which renders the phase diagram very 
complicated. Still unclear is the influence of the exter¬ 
nal field to the localized-to-critical transition. Moreover, 
the critical value of the spectral exponent was predicted 
by a recent mean-field analysis [22j to s* = 1/2, that 
stands at variance to the aforementioned DMRG result. 
It thus remains a challenging task to map out precisely 
the localized-to-critical phase transition as represented in 
Fig. mb) in the absence of an external field. 

In the absence of bias and tunneling, the two-bath 
model exhibits a high level of symmetry, inclu ding the 
parity symmetry and rotational symmetry |2lL I22l IA fit . 
In the localized phase, spontaneous symmetry breaking 
takes place due to strong spin-bath coupling. Hence, a 
symmetry analysis may help distinguish the critical and 
localized phases. In addition, a novel quantum phase 
transition from a doubly degenerate “localized phase” to 


another doubly degenerate “delocalized phase” is uncov¬ 
ered with respect to the ratio of the coupling strengths 
a/p within the two baths [22|. The transition is in¬ 
ferred to be of the first order, and the transition point 
a/P = 1 is determined when the spectral exponents of the 
two baths are identical. Since the system at the transi¬ 
tion point corresponds to the XZ-symmetric spin-boson 
model, the critical properties of the ground state, i.e., 
the spin polarization m(a, s) and generalized suscepti¬ 
bility y(a, s) at s = s and a = /3, can also be used to 
distinguish the localized and critical phases. 

The purpose of this paper is to investigate various 
ground-state phases in the extended spin-boson model 
involving two identical, independent baths, and deter¬ 
mine the critical value of the spectral exponent s* sepa¬ 
rating the localized and critical phases in the weak cou¬ 
pling limit of a —y 0. Via the variational approach, the 
DMRG approach, and the exact diagonalization method, 
we conduct a comprehensive study on the ground state 
properties of the two-bath spin-boson model with zero 
bias and tunneling for the baths described by single mode 
as well as continuous spectral density function. In this 
work, rotational and parity symmetry breaking is found 
to occur along the localized-to-critical phase transition, 
and the phase boundary is obtained with s* = 0.49(1), 
consistent with the mean-field predictions. 

The rest of the paper is organized as follows. In Sec. II, 
the two-bath spin-boson model and its symmetry prop¬ 
erties are described, along with an introduction to the 
variational approach. Sec. Ill and IV present the numer¬ 
ical results for the localized and critical phases in the 
two-bath spin-boson model involving a spin coupled to 
two oscillators or two baths described by a continuous 
spectral density function, respectively. Conclusions are 
drawn in the final Sec. V. 


II. METHODOLOGY 
A. Model 

The two-bath spin-boson model can be described by 
the Hamiltonian below 

e A 

H = -a z - —a x + 

+ -y i + fyi) 

l 

+ ~o~ Y, 2 + fy,2)? (1) 

l 

where e and A is the spin bias and tunneling constant, 
respectively, i = 1,2 is the index of the baths, and A; 
((pi) is the diagonal (off-diagonal) coupling strength. In 
order to investigate the quantum phase transition be¬ 
tween the critical and localized phases, we focus on the 
case of e = A = 0 as mentioned in the ‘Introduction’. A 
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logarithmic discretization procedure is adopted by divid¬ 
ing the phonon frequency domain [0, icJ into M intervals 
w c [A- z , (Z = 1,2,..., M) (37, 381- The coupling 

strengthes un and A/ (or pi ) in Eq. £[]) can then be cal¬ 
culated as 

/•A l ui c 

Xf= dtJ(t ), 

J A~ I ~ 1 u c 

uji = A; -2 / dtJ(t)t , (2) 

For convenience, the frequency cut off w c = 1 and the dis¬ 
cretization factor A = 2 are set throughout this paper. 
It should be noted that infinite bath modes are consid¬ 
ered via the integration of the continuous spectral density 
J(w), although the number of effective bath modes M is 
finite. 

Since various values of ( a x ) and (cr~) are possible due to 
the extended symmetry of the two-bath model with zero 
bias and tunneling, the spin polarization is introduced as 

m = yj (a x ) 2 + ( a y ) 2 + (cr 2 ) 2 (3) 

Due to Hamiltonian invariance under the transforma¬ 
tion <j y —> —<Jy, the y component of the spin polariza¬ 
tion is (a y ) = 0. Hence, m can be simplified to be 

Wx} 2 + Wz) 2 - According to Ref. J35| , the critical phase 
is characterized by (cq) =0 (i = x,y 7 z), which results in 
m = 0 . 

In addition, the phonon population P p h{x, z) is intro¬ 
duced to depict the ground state of the two-bath model. 
Assuming that the wave function of the ground state can 
be written as 

l^s) = l+)l^+)ph + 1“ )l'0-)ph, (4) 

where |'0+)ph an( 4 |'0-)ph are the phonon parts of the 
wave function corresponding to the spin up and down 
states, respectively, which can be expanded in a series of 
Fock states or coherent states. The phonon population 
P p h(x, z) is thus defined as 

Pph(x,z) {x, z) | ffg (x, z)) S pin 

= \^+{x,z)\ 2 + \ip-(x,z)\ 2 (5) 

where (• • • ) S pm represents the trace of the spin freedom in 
the wave function, x and z are coordinates in X — Z plane 
corresponding to the off-diagonal and diagonal coupling 
baths, and p±(x, z) = (f f l'0±)ph is the phonon-component 
of wave function in the two-dimensional coordinate rep¬ 
resentation r = ( x , z). 


B. Ground state symmetry 

The model studied in this paper exhibits a high level of 
symmetry due to zero bias and tunneling (e = 0, A = 0). 
A group theory analysis [2l|, [22j] shows that the ground 


state is always doubly degenerate. We introduce four 
parity symmetry operators including I = id as the unit 
operator, 

V x = a x e^^ b t^'\ ( 6 ) 

V z = a z e in ^‘ b l* bl ’ 2 , 


and VxVz- The influence of the parity symmetry oper¬ 
ations to the ground state G is displayed in Fig. El a). 
Under the operation V x {V z ), the sign of the coordinate 
values corresponding to the displacements of phonons in 
the diagonal coupling bath (off-diagonal coupling bath) 
will be changed. The symmetry parameters ( x and £ z of 
the parity symmetry are defined as 

C* = (* g \Vx\*g), 

Cz = ^ g \P z \* g ). (7) 

When a ^ P, the results ( z = 1 and Cc = 0 (Cz = 0 and 
C x = 1) are obtained for the localized phase (delocalized 
phase) in the two-bath spin boson [22j. However, ( x and 
in the case o f a = ( 3 are still unclear. A vanishing 
value of C = vCzT+Cz usually indicates breakdown of 
the parity symmetry. In contrast, one has £ = 1 when 
the ground state has perfect parity symmetry along the 
X or Z direction. If 0 < ( < 1, the ground state exhibits 
partial parity symmetry which may be induced by the 
numerical errors or the finite number of the degrees of 
freedom. 

In the XZ-symmetric spin-boson model with s = s and 
a = ft, the system may exhibit a rotational symmetry, 
since the Hamiltonian is invariant when one simultane¬ 
ously rotates the spin and the two baths in the X-Z plane 
by an arbitrary angle 6. According to the Abelian 1/(1) 
symmetry of the two-bath model proposed in Ref. 1361 ]. 
the rotational symmetry operator T(0) is introduced as 

T{0) = exp {-WS), (8) 

where S is the generator of the t/(l) symmetry defined 
as 

1 M 

£= 2 a v + { bl ^ lb U ~ b U bl ’ 2 ) ■ ( 9 ) 

;=i 

Figure. [U(b) shows the influence of the rotational symme¬ 
try operator on the ground state G, where the coordinate 
values of the ground state f x ,f z in the X and Z direc¬ 
tions are proportional to the displacement coefficients of 
the phonons in the off-diagonal coupling and diagonal 
coupling baths, respectively. 

The symmetry parameters y( 9 ) and 7 P h($) are intro¬ 
duced to quantitatively measure the rotational symmetry 
as 


7(0) 

7 P h(0) 


(* 9 |t(6>)|* 9 >, 

(^|f ph (0)|^) 


<*» I ^P 




( 10 ) 


!*«>■ 
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Figure 2: (Color on-line) The schematic of the influence of the 
parity symmetry operators V x ,Vz and rotational symmetry 
operator T(6 ) to the ground state G are displayed in (a) and 
(b), respectively. For convenience, we use the polygon shape 
of the ground states to emphasize the influence of symme¬ 
try operations. VxG,V z G and T(9)G are new ground states 
under these symmetry operations, and f x , f z represent the co¬ 
ordinate values of the ground state in the X and Z directions, 
respectively. 


It should be noted that only the rotation of the two baths 
is considered in the definition of 7 p h(0). The system 
energy E(9) = (4' g \T^(6)HT(9)\A> g ) is expected to be 
independent of the rotational angle 9 in the two-bath 
spin-boson model involving two identical baths. If the 
ground state has perfect rotational symmetry, one has 
7ph($) = 1 for the whole regime of rotational angle 9. 
In contrast, 7 P h($) ~ 5(9) is obtained in the absent of 
the rotational symmetry, where 5(6) is a delta function. 
When the ground state has partial rotational symmetry, 
7ph (9) decays with the rotational angle 9. Furthermore, 
we investigate the 0-dependent behavior of the parity- 
symmetry parameters (,(&) defined as 

Cr(0) = (* g \Tl h (9)V x f ph (9)\* g ) 

Cz(9) = {* g \f$ b {0)V t f ph (0)\* g ) 

m = VqW+ZW, (id 

where T p h(0) defined in Eq. (flOl) is the phonon part of the 
rotational operator, and T p h(0)|4' 9 ) is one of the degener¬ 
ate ground states obtained by rotating the ground state 
|4/ 9 ). In the rest of the paper, both symmetry param¬ 


eters, 7 P h($) and ((9) of the rotational and parity sym¬ 
metries, respectively, will be comprehensively studied, as 
they are useful and sensitive to detect the spontaneous 
symmetry breaking in the localized-to-critical phase tran¬ 
sition. 


C. Variational method 


A systematic coherent-state expansion of the ground 
state wave function, termed as the “multi-Di ansatz,” is 
introduced as the variational trial ansatz [22, |39|, j40| ■ It 
can be written as 


l*> 


N 

|+)y> ra exp 

n= 1 
N 

|-> ^2 Bn exp 

n— 1 




|0) P h 


H 


.c.) 


|0) p h(12) 


where H.c. denotes Hermitian conjugate, |+) (|—)) stands 
for the spin up (down) state, |0) p h is the vacuum state of 
the phonon bath, and M and N represent the numbers 
of the bath modes and coherent superposition states, re¬ 
spectively. In fact, Eq. m describes a superposition 
of the spin states |±) that are correlated with the bath 
modes with displacements f n j and g n ,i, where n and l 
represent the ranks of the coherent superposition state 
and effective bath mode, respectively. The displacements 
( fn,i,9n,i ) w ith 0 < l < M (M < l < 2 M) correspond 
to the phonons in the diagonal (off-diagonal) coupling 
bath. Using this trial wave function, the system en¬ 
ergy E can be calculated with the Hamiltonian expec¬ 
tation H = ('FlHj’]/) and the norm of the wave function 
D = (4'1'F) as E — H/D. The ground state is then ob¬ 
tained by minimizing the energy with respect to the vari¬ 
ational parameters A n , B n , f n j and g n j. The variational 
procedure entails N(AM + 2) self-consistent equations, 


dD = 

dxi dxi ’ 


(13) 


where Xi(i = 1, 2, • • • , ANM + 2N) denote the variational 
parameters. The “multi-Di” ansatz is much more sophis¬ 
ticated and contains more flexible variational parameters 
than the Silbey- Harris ansatz (4l| and Nazir’s ansatz [42f . 
where only 2 M + 1 and 4 M + 2 variational parameters 
are employed, respectively. For example, if N = 16 and 
M = 20, our new ansatz has 1312 variational param¬ 
eters, compared to 41 parameters in the Silbey-Harris 
ansatz and 82 parameters in Nazir’s ansatz. 

For each set of the coefficients (a, /?, s and s) in the 
continuous spectral densities J x (oj) and J 2 (w), more 
than 100 initial states are used in the iteration proce¬ 
dure with variational parameters (A n , B n ) uniformly dis¬ 
tributed within an interval [—1,1]. Displacement coef¬ 
ficients ( f n ,i,9n,i ) of the initial states obey the classi¬ 
cal displacements, i.e., f n j = —g n ,i ~ A// 2 u 7 for the 
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diagonal coupling bath and f n> i = —g n ,i ~ <^/ 2 w; for 
the off-diagonal coupling bath. In the single-mode case, 
fn,i , 9n,i , A n and B n are all initialized randomly. After 
preparing the initial state, a relaxation iteration tech¬ 
nique I Id. 11 ! is adopted, and a simulated annealing al¬ 
gorithm [ 22 j| is also employed to improve the energy min¬ 
imization procedure. The iterative procedure is carried 
out until the target precision of 1 x 10 -12 is reached. 

Theoretically, the number of coherent superposition 
states N —> oo is required for the completeness of the en¬ 
vironmental wave function in variational method. How¬ 
ever, large values of N pose significant challenges in car¬ 
rying out numerical simulations. To obtain reliable nu¬ 
merical results with large TV, an approach to improve the 
variational algorithm is undertaken based on the parity 
symmetry. Assuming |^ 9 ) is the ground state obtained 
by the variational method with TV coherent superposition 
states, an intermediate state I'l'int) can be generated via 
the parity symmetry operators I,V X ,V~ and V X V Z , 

|*int> = {Cil + C 2 V x + C 3 V Z + CiV x V z ) |tf 9 >, (14) 

where Ci(i = 1 , 2 ,3,4) is the weight coefficient. Accord¬ 
ing to the symmetry analysis [ 22 |, the symmetry oper¬ 
ator V x or V z can lead to the other branch of the dou¬ 
bly degenerate ground state with the same energy E g . 
Hence these four symmetry operators should be equally 
weighted with |Ci| = \C 2 \ = |C 3 | = |C 4 |. If Ci = 1, 
C 2 = ± 1,63 = ±1 and C 3 = ±1 can be derived. Simi¬ 
lar to the case of the delocalized Davydov D\ variational 
ansatz in the Holstein model [45}, the energy E { nt of the 
intermediate state is lower than Eg after considering 
the parity symmetry. Using these eight states as ini¬ 
tial states, one can obtain a new ground state | v I/ 9 ) by 
performing the variational procedure with 47V coherent 
superposition states, which yields a lower ground state 
energy Ef < E int < Eg. 

Due to numerical errors, however, the state T found 
by the variational algorithm corresponds only to the local 
minimum in energy in the vicinity of the ground state. 
To refine the variational results, the rotational symme¬ 
try should also be considered in the case of s = s and 
a = f3. Via the rotational operator T(6) acting onto the 
state | T), a subspace composed of a series of states with 
respect to the rotational angle 9 is built. Subsequently, 
the state with the minimum energy in this subspace is 
regarded as the ground state. Since the generater of the 
£7(1) symmetry S involves a hopping between the diago¬ 
nal and off-diagonal coupling baths, the displacement co¬ 
efficients in the two baths are identical after considering 
the rotational symmetry, consistent with the argument 
that the ground state is accompanied by a symmetric 
distribution of phonon numbers in the diagonal and off- 
diagonal coupling baths. 

The phonon population P p h(x, z) in Eq. (J5J can be cal¬ 
culated with the multi-Di variational ansatz in Eq. m 


as 

P ph (x, Z ) = J2 ^ AnF ^ x ' z ^' + ^ BnGn ^ x ' z ^\ ( 15 ) 

n =1 

where D = (Tgl^/g) is the norm of the wave function, A n 
and B n denote the weight coefficients of the n-th coher¬ 
ent superposition state coupled to the spin up and down 
states, and the phonon functions F n (x,z) = (r|'i/>+) p h = 
fn,x{ x )f™,z{ z ) and G n {x,z) = (r|^_) ph = gn, x {x)g n A z ) 
represent the phonon component of the wave function 
|'0±)ph in the two-dimensional coordinate representation 
f = (x,z). The function fn. x {x) denoting a coherent 
state in the off-diagonal coupling bath can then be de¬ 
duced, 


fn,x(x) = J|(x| /„,*) 

l 

(16) 

— ^ e~ ixiPl ^ 2 e ipix e~ u} ^ x ~ 

1 77 

- x ,) 2 /2 

where xi and pi are defined as 

Pi = 

(17) 

Xl - nr— ( fn,l + f n< l ) ■ 

V 2 loi 

(18) 


In the same way, the functions f n ,z(z), g n<x (x) and g ntZ (z) 
can also be calculated with the displacement coefficients 
fn,l and g n j in Eq. m as input. In the single-mode 
case, i.e., M = 1, the phonon function can be simplified 
as F n (x,z) = fn,x(x)fn, z (z ) = {z\f n ,i){ x \fn,2) where the 
subscripts 1 and 2 correspond to the diagonal and off- 
diagonal coupling oscillators, respectively. 


III. SINGLE MODE 

The ground state of the model involving two oscillators 
coupled diagonally and off-diagonally to a spin is inves¬ 
tigated in this section. The corresponding Hamiltonian 
can be written as 

-^single = w(b\bi + b\b 2 ) +—\(b\ + bi) (19) 

+ _ ^ _< ( , (^2 + ^ 2 ), ( 20 ) 

where A and <f> are diagonal and off-diagonal coupling 
constants, respectively. It is the simplest version of the 
two-bath model under current study. Furthermore, we 
focus on the case of two identical coupling constants A = 
<j) as it gives the Hamiltonian the rotational symmetry, 
which may provide some simple insights on the nature 
of the phase transition between the critical and localized 
phases. 
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A. Exact diagonalization 


The Hamiltonian is then described by 


In the exact diagonalization procedure, the phonon 
states |'0_|_) p h and |i/>_) p h corresponding to the spin up 
and down states, respectively, are expanded in a series of 
Fock states, 


N t r 

iv , +)ph = ^2 ( 2i ) 

kik2 


Ntr 

IV 1 —)ph = (22) 

k\k2 

where Ck 1} k 2 and dk Xl k 2 are the coefficients of the Fock 
state \k 1 k 2 ) for the two oscillators coupled diagonally 
and off-diagonally to the spin, and N tI = 100 is the 
bosonic truncation number defined as the cutoff value 
of the phonon occupation number. We have verified that 
this value of N tI is sufficiently large for the ground-state 
energy to converge. Solving the Schrodinger equation, 
one can obtain the wave function of the ground state 
|^/ s ) with a series of coefficients, Ck 1: k 2 and dk l7 k 2 , and 
the ground state energy E g . Thus, the expectation val¬ 
ues of a z and a x can be calculated as 

Ntr 

i&x) = ^ W fci k 2 d kl k 2 + d kl Ck-t k 2 1 

k\k 2 
Ntr 

( a z) = \ c kik 2 | — \dktk 2 | ■ (23) 

ki k 2 


The phonon population P p h(x, z) can also be obtained 
by 


Pph(x,z) = (24) 

Ntr 

y ( ^[Cfci,fe 2 Cfci,fc 2 (Xj ")] T [dk 1 ,k 2 Dk 1 ,k 2 {x, z)\ ^ , 
k\ k 2 


where Ck lt k 2 (x, z) and Dk lt k 2 (x, z) represent the Fock 
states \k 1 k 2 ) in the coordinate representation r = (x,z), 
corresponding to the spin up and down states, respec¬ 
tively. 


B. Analytical results 

The characteristics of the single-mode spin-boson 
model involving two oscillators can be investigated in¬ 
tuitively in coordinate representation with the transfor¬ 
mation x = (bi + /y/2 iw and z = (b 2 + /y/2ui. In 

the following discussion, we use x and z as the classi¬ 
cal counterparts of the corresponding operators x and z. 


H = H 0 + V 7 

H 0 = -iv 2 + i W V, (25) 

V = X'f-a, 

where we denote r = (x, z), r = \/x 2 + z 2 , a — (a x , <j z ), 
X' = Aa/w/ 2, and V 2 is the two-dimensional Laplace op¬ 
erator written in polar coordinates. This Hamiltonian 
describes a spin in a two-dimensional harmonic potential 
with spin-orbital coupling V. When A 3>w, the spatial 
motion of the particle is too slow compared to the degree 
of freedom of the spin. Therefore, it is justifiable to intro¬ 
duce the Born-Oppenheimer approximation. The spatial 
motion is thus treated classically, and the spin-part is 
solved by 


V\r]±) = e±\rj ± ) 


with the adiabatic eigenstates 


\v+) 

\v~) 


"cos(f - 
sin(f - 

sin (f _ 

- cos (f 


f) 

I) J’ 

-!) 
-§)J 


(26) 


(27) 


and the corresponding eigenvalues e± = ±AV. Here we 
define tand = z/x. 

Then, the wave function of the system can be assumed, 


\V) = <P+(f) \v+) + <P-(r) \V~) ■ (28) 


Using the adiabatic eigenstate as the basis, the equations 
for the spatial part of the wave function are obtained with 
the stationary Schrodinger equation H |T) = E |\H), 


. V 2 1 
(_ ^“ + 8 ^ + 
= 0 , 


2 2 
urr* 


+ X'ra z +iV n - a a v 


E)$(r) 

(29) 


where E = diag {E + ,E_}, #(f) = <^_(r)) T , <x z 

and a y are the Pauli matrixes, and the non-adiabatic 
terms are given by the operator 


yn.a. 


J_8_ 

2r 2 86 ’ 


(30) 


which are only in the angular direction, and can be ne¬ 
glected in further analysis. Following the standard pro¬ 
cedure of variables’ separation method, the solution has 
the following form 


<P±{r)='^2 Cme me R± (r, m ), (31) 

m 


wherein the radial function R± (r, m) is determined by 
the equation 


' 1 8 

( r d \ 

2 r dr 

iy 8r) 



R± (r, to) 


E±R± (r, to) 
(32) 
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Figure 3: (Color on-line) The wave function of the ground 
state in a strong coupling case of A = (j> = 10 and w = 1 is 
displayed in two-dimensional coordinate representation (x,z). 
The x- and ^-coordinate correspond to off-diagonal and diag¬ 
onal coupling baths, respectively, and the colour represents 
the phonon population Pph(x,z). In (a) and (c), the numbers 
of coherent superposition states N = 8 and 32 are adopted, 
respectively, and an intermediate state defined in Eq. m is 
shown in (b). 
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Figure 4: (Color on-line) The convergence test of the ground 
state at A = <j> = 10 and w = 1 is displayed for various num¬ 
bers of the coherent superpositions N = 4,6, 8, 12, 32 and 96 
in (a)-(f), respectively. The x- and ^-coordinate correspond 
to off-diagonal and diagonal coupling baths, respectively, and 
the colour represents the phonon population P p h{x, z). 


with the effective potential 


y eff( r ) = Y r2±Xr+ (g 



(33) 


The effective potential contains a harmonic potential, a 
linear potential and a centrifugal potential, and the an¬ 
gular quantum number is half-integer (m = ±|, ±|...) 
due to the contribution of the spin 1/2 part. Thus, the 
ground state that corresponds to m = ±1/2 is doubly 
degenerate. In the strong coupling regime (A to), we 
can neglect the centrifugal potential, leading to the ex¬ 
pectation of a ring-shaped ground state with the radius 
R oc X/uj 2 . 


C. Numerical results 

We first investigate the ground state of the two-bath 
model in the case of w = 1 and A = <p = 10 via the 
variational method with N = 8 and M = 1. According 
to the aforementioned theoretical arguments, the wave 
function of the ground state is expected to be ring shaped 
in the two-dimensional ( x , z) coordinate representation. 
However, variational results depict only a quarter of the 
ring as shown in Fig. [U)a), where the colour represents 
the value of the phonon population P p h(x, z) defined in 
Eq. (fl5l) , and x- and z-coordinates correspond to the off- 


diagonal and diagonal coupling baths, respectively. Us¬ 
ing the parity symmetry operators onto the ground state, 
an intermediate state defined in Eq. m is obtained and 
shown in Fig. [U)b) . The energy of this intermediate state, 
Pint = —25.48911811, is found to be slightly below the 
ground state energy P^ =s = —25.48058088. Taking this 
intermediate state as an initial state, one can seek the 
ground state via the variational method with N = 32. 
Figure He) shows the phonon population P p h(x, z) de¬ 
creasing smoothly with the z-coordinate, quite different 
from that of the intermediate state. It indicates that 
the parity symmetry in the z direction is broken, re¬ 
sulting in doubly degenerate ground states with differ¬ 
ent values of ( a z } but the same ground state energy 
Eg~ 32 = —25.49741979 that is much lower than both 
Eg =8 and P in t- 

Further, the convergence of ground state with respect 
to the number of the coherent superposition states N 
warrants a careful examination. As shown in Fig. U the 
phonon population P p h(x, z) gradually starts to resemble 
a ring-like shape as N is increased from 4 to 6,8,12, 32 
and 96. In fact, the X-Z symmetric spin-boson model ex¬ 
hibits continuous degeneracy by the projecting T(0)|'F g }, 
where T(6) is the rotational symmetry operator defined 
in Eqs. (0) and and |T g ) is one branch of the ground 
state. Moreover, the ground state energy E^ monoton- 
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Figure 5: (Color on-line) The wave function of the ground 
state in a case of A = 4> = 2 and w = 1 is displayed in (a)-(d) 
for various numbers of the coherent superposition states N = 
4, 6,16 and 24. The x- and ^-coordinate correspond to off- 
diagonal and diagonal coupling baths, respectively, and the 
colour reflects the value of the phonon population P p h(x,z). 


ically decreases with N and is convergent to an asymp¬ 
totic value Eg~ 96 = —25.497421539, consistent with the 
exact diagonalization result E g = —25.497421544 with 
the phonon truncation number N tI = 100. 

The ground state of the two-bath model in the case of 
weaker coupling, A = </> = 2 and u = 1, is investigated 
next. The phonon distribution P p h(x,z) for N = 4,6,16 
and 24 is displayed in Fig. El)a)-(d). Different from the 


0.31 20 


I 


(a) ■ (b) 



Figure 6: (Color on-line) The wave function of the ground 
state obtained by variational method is displayed in (a)-(d) for 
the coupling strengths A = (j> = 0.1, 2,10 and 20, respectively. 
The phonon frequency u = 1 is set for both two baths. The 
x- and z-coordinate correspond to off-diagonal and diagonal 
coupling baths, respectively, and the colour represents the 
phonon population P p h(x,z). 


Figure 7: (Color on-line) The wave function of the ground 
state obtained by exact diagonalization is displayed in (a) and 
(c) for A = <j> = 10 and A = <j> = 2, respectively. Correspond¬ 
ingly, the difference between the exact diagonalization and 
variational results are displayed in (b) and (d). The x- and z- 
coordinate correspond to off-diagonal and diagonal coupling 
baths, respectively, and the colour represents the phonon pop¬ 
ulation P p h(a;, z). 


results shown in Fig. |4] the shape of the ground state 
remains nearly unchanged, indicating that a small value 
of N is sufficient to obtain a reliable numerical result. 
The ground state energy £)^ =24 = —1.368929967 is again 
in an excellent agreement with the exact diagonalization 
result E g = —1.368929970. 

Finally, the ground states of two-bath model with 
A = 4> = 0.1 and 20 are also plotted in Fig. [6])a) and 
(d), respectively, to facilitate comparison with those for 
A = <j> = 2 in Fig. [Hlb) and A = <j> = 10 in Fig. (HDc). In 
the weak coupling regime A -C u, a ground state with a 
clear rotational symmetry is found, while it collapses in 
a corner of the X-Z plane in the strong coupling regime 
A to. It supports our conjecture that the rotational 
symmetry breaks when the coupling strength exceeds a 
certain value A c , similar to the picture of the phase tran¬ 
sition in the classical XY model. Since the radius of the 
circle in Fig. [6]) a) is quite small, any slight shift of the cen¬ 
ter from the coordinate origin (0, 0) will induce a sharp 
jump in the spin polarization from m = 0 to m k ±1. It 
indicates that the spin polarization m is unstable in the 
weak coupling regime, corresponding to the free phase. 
That the ground state in Fig. Etc) shows a crescent profile 
rather than a complete ring may be reasoned from pre¬ 
vious analytical results. The ground state here must be 
doubly degenerate, and our numerical calculations yield 
only one branch of the ground state. Upon combining 
both the degenerate sates, once can readily obtain the 
complete ring shape of the ground state. 
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Table I: The ground state energy E g and spin polarization 
m obtained by the variational method (VM) and exact diag- 
onalization (ED) are listed for three different cases with the 
diagonal and off-diagonal coupling strengths A = (j> = 0.1,2 
and 10, respectively. The phonon frequency uj = 1 is set, and 
N and Ntv denote the coherent-superposition number in the 
variational method and the bosonic truncated number in the 
exact diagonalization, respectively. 




II 

-O- 

II 

o 

i— 1 

CN 

II 

II 

A = (/>= 10 

VM 

N 

8 

24 

96 


E s 

—4.987582654E—3 

-1.368929967 

-25.497421539 


m 

0.995049457657 

0.568606150 

0.500053226 

ED 

N tI 

100 

100 

100 


E s 

—4.987582654E—3 

-1.368929970 

-25.497421544 


m 

0.995049457657 

0.568606143 

0.500053212 


D. Discussion 

The ground sates obtained by the exact diagonaliza¬ 
tion method are shown in Figs.[2 a ) and [2c) for the two 
cases of A = </> = 10 and A = <j> = 2, respectively. Results 
for both of these cases seem to be identical to those with 
the variational results shown in Figs. H]T) and (2d). To 
further verify the consistency, the difference between ex¬ 
act diagonalization and variational results are displayed 
in Figs. (2b) and (d). The resulting difference is two or¬ 
ders or magnitude smaller than the phonon population 
Pph(x,z), implying that the wave function obtained by 
the two methods are nearly the same, thereby lending 
support to the superior accuracy of our variational re¬ 
sults. 

The ground state energy E g and spin polarization m 
obtained by the variational method and the exact diag¬ 
onalization approach are summarized in Table. Q] for the 
coupling strengths A = <f> = 0.1,2 and 10. In all three 
cases, results from both the methods agree to each other 
for more than 9 significant digits of E g and m. More¬ 
over, the radii of the rings R = 5,25 and 50 in the cases 
of A = (j) = 2,10 and 20 calculated from Figs. 032), [4))f) 
and (2d), respectively, are found to be consistent with 
the theoretical prediction R= cX/oj 2 with the coefficient 
c = 2.5. This excellent reproduction of results again 
points to the superiority of the variational method and 
to the robustness of the ground state obtained by numer¬ 
ical calculations. 

The symmetry of the ground state in the single-mode 
case is also studied via the symmetry parameters ( = 
a/C-c + Cz °f the parity symmetry and y p h of the rota¬ 
tional symmetry. Though the phase transitions may be 
reduced to the ground-state level crossings due to the 
finite number of degrees of freedom, the symmetry prop¬ 
erties of the ground states in the localized (A = (j) 1), 

critical, and free phases (A = <f> <C 1) are unchanged. 
Fig. [8] shows ((0) for A = cjj = 0.1,2,10 and 20 when 


s = A = 0 and uj = 1. Interestingly, it is found that ( = 1 
regardless of 9 in the weak and intermediate coupling 
regimes, pointing to the parity symmetry in the critical 
and free phases. In the localized phase, however, narrow 
peaks of ((6) are found for a strong coupling strength 
A = 20. It indicates that the ground state is localized 
without the parity symmetry. In Fig. (9[ the symmetry 
parameter 7 P h($) of the rotational symmetry is also plot¬ 
ted, which shows an abrupt decay to zero in the strong 
coupling regime (A = <j> = 20) but remains equal to one 
in the weak coupling regime (A = <jj = 0.1). In the in¬ 
termediate regime, 7 P h($) is found to decrease gradually. 
These numerical results further support our contention 
that the rotational symmetry breaks only when the cou¬ 
pling is strong. 

IV. CONTINUOUS SPECTRAL DENSITIES 
A. The case with a = /? 

In this subsection, we study the ground state proper¬ 
ties of the two-bath model involving the baths described 
by a continuous spectral density function J(ui) via the 
variational approach. Infinite bath modes are considered 
in the variational calculations, although the number of 
the effective modes M is finite in the logarithmic dis¬ 
cretization procedure. For convenience, we first examine 
the case involving two identical baths, i.e, s = s and 
a = fi. 

The spin polarization m defined in Eq. m is displayed 
in Fig.[[n]as a function of the coupling strength a for var¬ 
ious values of the spectral exponent s = 0.4, 0.5,0.6, 0.7 



Figure 8: (Color on-line) The parity symmetry parameter 
m = Va + ££ is displayed as a function of the rotation 
angle 8/n for coupling strengths A = (j>. The spin bias e = 
0, tunneling constant A = 0, and frequency uj = 1 are set. 
The dash lines represent the fitting with the trigonometric 
functions. 
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Figure 9: (Color on-line) The rotation symmetry parameter 
7 P h(0) is displayed as a function of the rotation angle 6/n for 
various coupling strengths A = 4> = 0.1, 2,10 and 20. The spin 
bias e = 0, tunneling constant A = 0, and frequency w = 1 
are set. 



Figure 10: (Color online) The spin polarization m(a) at var¬ 
ious values of the spectral exponent s is plotted as a function 
of the coupling strength a in the case with a = /? and s = s. 
The number of coherent superposition states N = 4 and effec¬ 
tive bath modes M = 20 are used in variational calculations. 
The downward arrows indicate the transition points a c . 


and 0.8 in the case of s = s and a = /3. For a > a c , 
an increase of the spin polarization m(a) is found for 
all of s, corresponding to the localized phase shown in 
Fig. mb). However, a non-zero spin polarization is found 
in the critical phase with a < a c , quite different from the 
prediction of m = 0 by an earlier study [35j . In addition, 
the localized-to-critical transition point a c marked by the 
downward arrows is shifted visibly with an increase in s 
except for s = 0.4 for which no phase transition occurs. 
It indicates the critical value of the spectral exponent 
is s* ~ 0.5, consistent with the prediction s* = 1/2 


Figure 11: The ground state energy Eg is displayed as a 
function of the coherent-states number A in a case of the 
critical phase with s = s = 0.8 and a = 0 = 0.2. The 
number of the effective bath modes M = 20 is used in vari¬ 
ational calculations. The dash-dotted line represents the fit¬ 
ting Eg = aN~ b + £ g (oo). In the inset, the spin polarization 
m(N) is also shown on log-log scale, and the dashed line in¬ 
dicates a power law fit. 


of the mean-field analysis [22| . but much smaller than 
s * = 0.75(1) m. 

The convergence of the variational results with re¬ 
spect to N and M is carefully tested. Fig. [TT] shows 
the ground-state energy Eg as a function of N in the 
critical phase with s = s = 0.8 and a = /3 = 0.2. A 
power law decay of the ground state energy with the 


- 0.1918 



- 0.19185 


- 0.1919 


- 0.19195 



0.04 0.06 0.08 0.1 

1/M 


Figure 12: The ground state energy Eg and the spin polar¬ 
ization m(M) are displayed as a function of the number of the 
effective bath modes M at N = 4 in a case of the critical phase 
case with s = s = 0.8 and a = /? = 0.2. The dash-dotted and 
dashed lines represent the fitting y(M) = aM~ b + y( oo). 
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form = aN~ b + E g (oo) is found via numerical fitting, 
which yields the asymptotic value E g ( oo) = —0.193572. 
In the inset, the spin polarization m(N) is also displayed 
as a function of N on a log-log scale. A perfect power- 
law behavior of m(N) is obtained with the slope 0.174(2) 
and the asymptotic value m( oo) = 0. It suggests that 
the non-zero value of the spin polarization in the critical 
phase originates in the effects of the finite N. Further¬ 
more, the spin polarization m(N) in the critical phase is 
not convergent even for N = 96, unlike the case of the 
localized phase where a small value of N is sufficient to 
obtain reliable results. In a similar manner, the influence 
of M to the ground state energy E ^ 1 a nd spin polariza¬ 
tion m(M) is also depicted in Fig. 1121 Both quantities 
are found to reach asymptotic values when 1/M < 0.05, 
indicating the sufficiency of M = 20. Therefore, in the 
following discussion on the variational results the number 
of coherent superposition states and the bath modes are 
set to TV = 16 and M = 20, unless specified otherwise. 

Table. El presents a comparison between the numeri¬ 
cal results obtained bv the variational method and the 
DMRG approach [37l I46L l47j|. To ensure reliable results, 
the phonon number used in DMRG algorithm is d p = 50, 
much larger than d p = 30 used in the previous work [35j . 
The length of Wilson chain is set to L = 50 and the cutoff 
dimension of the matrix is D c = 60. In the three cases 
of (s = s = 0.4, a = /3 = 0.1), (s = s = 0.6, a = /? = 0.1) 
and (s = s = 0.8, a = /? = 0.2), only a slight difference 
of the ground state energy, i.e., A E/E g < 0.1%, is found 
between the variational and DMRG results, further re¬ 
inforcing the superiority of our variational results. Since 
the ground state of the critical phase is unstable I35| . 
the spin polarization obtained by the variational method 
is larger than that by the DMRG, as shown in the last 
two columns in Table. El In the localized phase, how- 


Table II: The ground state energy E g and spin polarization 
m obtained by the variational method (VM) and DMRG are 
listed for three different cases with (s = s = 0.4, a = (5 = 0.1), 
(s = s = 0.6, a = P = 0.1) and (s = s = 0.8, a = j3 = 0.2) 
in the localized phase (first case) and critical phase (last two 
cases). Three numbers of the coherent superposition states 
N = 16, 64 and 96 are used in variational calculations, which 
are sufficiently large in each case. d p = 50 represents the 
phonon number allocated on each site on the Wilson chain in 
the DMRG algorithm. 




s = s = 0.4 

a = /3 = 0.1 

s = s = 0.6 

a = 13 = 0.1 

s = s = 0.8 

a = 13 = 0.2 

VM 

N 

16 

64 

96 


E s 

-0.16759 

-0.12917 

-0.19356 


m 

0.77448 

0.53372 

0.40171 

DMRG 

dp 

50 

50 

50 


E s 

-0.16771 

-0.12923 

-0.19357 


m 

0.75496 

0.42665 

0.29635 



Figure 13: (Color on-line) The spin polarization m and its 
x and z components (a x ) and (a z ) are plotted as a function 
of the rotational angle 9 for the states T(0)|'h g ) in the case 
of s = s = 0.8 and a = /.3 = 0.2. The dash lines represent 
the fitting with the trigonometric functions. In the inset, the 
shift A E = E{6) — E g from the ground state energy is shown. 



Figure 14: (Color on-line) The symmetry parameter of the 
parity symmetry (( 0 ) — ^/( x ( 9) 2 + ( z ( 9) 2 is displayed as a 
function of 9/ n for the two cases of s = s = 0.4 (bottom) and 
s = s = 0.85 (top) when the coupling strengths are a = /? = 
0.02. The width of the peak A 9 is defined as the size of the 
parity-symmetry regime with £ > 0. 

ever, a small value of N = 16 is sufficient to obtain the 
variational result of the spin polarization m = 0.77448, 
comparable with the DMRG result m = 0.75496. 

Figure. fl3l shows the spin polarization m(0) and its x 
and 0 components (a x ) and (tr z ) for the states T(0)|'F s ), 
where T{ff) is the rotational symmetry operator defined 
in Eq. © and | H /q ) is the ground state obtained by 
the variational method. As the rotational angle 9 in- 
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Figure 15: (Color on-line) The parity index, P(a, s) = 

2Ad/n, is displayed in (a) as a function of s for various val¬ 
ues of a and in (b) as a function of a for various values of 
s. The transition point separating the localized phase from 
the critical phase is located at the position where the parity 
index P = 1 reaches. 


creases, the values of (a x ) and ( a z ) oscillate between 
—0.5 and 0.5, while the corresponding spin polarization 
m remains almost unchanged. The obtained curves can 
be fitted with trigonometric functions, ( a x ) = —msin(0) 
and (a z ) = mcos(0) with m = 0.53091. It indicates that 
neither ( a x } nor (cr 2 ) is a good candidate to character¬ 
ize the localized-to-critical phase transition, even though 
they were employed in Ref. [35}. In the inset, the shift 
A E(9) = E(9) — E g is plotted. The sufficiently small 
value of A E k 2x 10 -15 shows that there are continuous 
degenerate ground states which have the same energy E g , 
consistent with the prediction from the rotational sym¬ 
metry analysis. 

Apart from the spin polarization, we have also probed 
the symmetry properties of the ground state. The 
symmetry parameter of the parity symmetry, £(0) = 
VCx(0) 2 +( z (6) 2 where ( x and ( z are defined in Eq. 0, 
is displayed in Fig. [14] for the two cases of s = s = 0.4 
and 0.85 at a sufficiently small coupling strength a = 
f3 = 0.02. In the localized phase with s = 0.4 < s*, 
sharp peaks of £(0) are found at 9/tt ss n/2 {n = 1,2,3 



Figure 16: (Color on-line) The symmetry parameter of the 
rotation symmetry 7 P h($) is displayed against the rotational 
angle 8/n for the two cases of s = s = 0.4, a = /3 = 0.02 and 
s = s = 0.85, a = P = 0.02, corresponding to the critical and 
localized phases, respectively. 


and 4) with a small peak width A 9 defined as the size 
of the parity-symmetry regime ( > 0. It suggests that 
the ground state is localized in a corner of the X-Z plane. 
On the other hand, £(9) is always greater than zero in 
the critical phase with s = 0.85 > s*, indicating that the 
parity symmetry covers the whole subspace. Hence, the 
parity index P = 2A9/tt reflecting the localization of the 
ground state can be used to quantitatively distinguish 
the localized and critical phases. 

Shown in Fig. HU a) is the parity index P(a, s ) versus 
the spectral exponent s for various values of a, in the case 
with two identical baths, i.e., s — s and a = / 3 . Without 
any loss of generality, it is assumed that only the ground 
state with the parity index P = 1 belongs to the critical 
phase, otherwise it belongs to the localized phase. Ac¬ 
cording to this criteria, the transition point s c between 
the localized and critical phases is calculated for various 
values of a. With an increase in a, s c increases monoton- 
ically, in agreement with the trend shown in Fig. [ljlo) . 
Moreover, the parity index P(a, s) is also displayed in 
Fig- EEM b) for various values of s. The transition point 
a c can be measured as a function of s in a similar manner 
and subsequently the phase boundary in the X-Z plane 
can be obtained. 

The rotational symmetry is studied next in the criti¬ 
cal and localized phases, using the two typical cases of 
s = s = 0.4, a = j3 = 0.02 and s = s = 0.85, a = /3 = 
0.02, respectively. The symmetry parameter 7 P h(0) is 
displayed in Fig. [16] by rotating the two baths in the X-Z 
plane through an angle 9. In the localized phase, 7 p h($) 
quickly depletes to zero, different from that in the criti¬ 
cal phase, where it gradually decays to a nonzero value. 
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Figure 17: (Color on-line) The symmetry parameter of the 
rotational symmetry 7 P h(a, s) at the angle 6 = 1/27r is dis¬ 
played in (a) as a function of s for various values of a and in 
(b) as a function of a for various value of s. The transition 
point is located at the position separating the localized phase 
with 7 P h = 0 from the critical phase with 7 P h > 0. 


Compared to the results of the single-mode case in Fig. 0 
one can find that the localized and critical phases corre¬ 
spond to the strong and intermediate coupling regimes, 
respectively. In general, the strong coupling regime has a 
large coupling strength a and a small spectral exponent 
s, while opposite trends ensue in the intermediate cou¬ 
pling regime. For the spectral exponent s > 1, however, 
the system always resides in the weak coupling regime, 
corresponding to the free phase. 

Similar to the case of the parity symmetry, the sym¬ 
metry parameter of the rotational symmetry 7 p h(a,s) 
can also be used to distinguish the localized and critical 
phases. Without loss of generality, we set a special rota¬ 
tional angle 6 — 7r/2 where the phonons in the diagonal 
and off-diagonal coupling baths are interchanged. Thus, 
the value of the symmetry parameter is expected to be 
7 ph = 0 in the localized phase and 7 P h > 0 in the critical 
phase. Fig. fTTT a) shows 7 P h(s) as a function of s for vari¬ 
ous values of a. The transition point s c is then located at 
the position separating the localized phase from the crit¬ 
ical phase. It monotonously increases with a, consistent 
with the results in Fig. [151 a). Moreover, 7 P h(<a) is also 



Figure 18: The rotation symmetry parameter 7 P h(a, s) and 
the parity index P(a, s) for TV = 2,12,16 and 20 are displayed 
as a function of s in a weak coupling case of a = /3 = 0.01. 
The dash line indicates the transition point s c = 0.495(5). 

plotted in Fig. fTTI b) as a function of a for various val¬ 
ues of s. The transition boundary a c (s) separating the 
localized and critical phases can then be appropriately 
calculated. 

To accurately estimate the critical value of the spec¬ 
tral exponent s* shown in Fig. |ljb), the case of very 
weak coupling strength a = /3 = 0.01is used to inves¬ 
tigate 7 P h(a, s) and P(a,s). As shown in Fig. [TSl tran¬ 
sition behavior of P(a,s ) and 7 P h(a, s) is displayed for 
TV = 2,12,16 and 20. The near overlap of the two curves 
of TV = 16 and 20 suggests that TV = 16 is sufficiently 



Figure 19: Displayed as a function of the superposition num¬ 
ber TV is the transition points determined by the variational 
calculations at weak coupling a = /3 = 0.01 (lower curves) 
and 0.02 (upper curves). The circles and triangles correspond 
to the results from the curves of 7 P h(a, s) and P(a, s), respec¬ 
tively. The dash lines are the linear fitting for the extrapola¬ 
tion of s c with 1/TV — > 0. 
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large to accurately describe the phase transition. With 
the increase of N, the transition point on the P(a, s) and 
7 P h(a, s ) lines is found to decrease monotonically tending 
to an asymptotic value of s c = 0.495(5) as marked by the 
dashed line. It points to the critical value of s* = 0.49(1) 
in the weak coupling limit of a 0, in perfect agree¬ 
ment with the mean-field prediction of 1/2, but stands 
at variance with the value of 0.75(1) by Guo et al. )35j] . 
This discrepancy may be attributed to the fact that the 
numerically unstable critical phase is beyond the reach of 
the DMRG algorithm of Guo et al., and the external field 
holds great sway over the phase-transition properties of 
the two-bath model. 

In order to get a good estimate of s c (N), transition 
points calculated from the variational approach are in 
Fig. [T9] as a function of 1/N where N is the number 
of superpositions. Two values of coupling strengthes 
a = /3 = 0.01 and 0.02 are used, and we also set 
e = A = 0, and s = s. As 1/N decreases, the differ¬ 
ence between 7 P h(<a, s) and P(a, s ) gradually disappears 
for both values of a. Using linear fitting of s c (l/N), 
the asymptotic values of s c = 0.493(6) and 0.541(7) for 
the two cases are obtained by extrapolation to infinite N, 
which is consistent with the N = 16 results s c = 0.500(5) 
and 0.55(1) within the error bars, further supporting our 
claim that the number of superpositions N = 16 is suf¬ 
ficient to obtain reliable results. The deviation of the 
critical point from the Guo’s result s* = 0.75(1) is not 
induced by the effect of the finite value of N. The phase 
diagram of the extended spin-boson model is displayed 
in Fig. (TU] in the case of s = s and a = (3. The solid tri¬ 
angles and stars represent the phase boundary obtained 
from the parity index P(a, s ) and symmetry parameter 
of the rotational symmetry 7 P h(a, s), respectively. The 
error bars in this phase diagram are estimated via the dif¬ 
ference of the transition points measured with the fixed 
a and s, respectively. 

Having studied the symmetry properties in detail, we 
now turn our attention to the wave function of the ground 
state for the two-bath model involving the continuous 
spectral density to understand the nature of the local¬ 
ized and critical phases. To serve this purpose, we chose 
the case of s = 0.5, a = 0.9 in the localized phase and 
the case of s = 0.85, a = 0.02 in the critical phase as 
examples. Fig. [TT] shows the displacement coefficients 
fn,i and g n g defined in Eq. (fl2l) in the localized phase 
as a function of the phonon frequency wj. Two differ¬ 
ent behaviors of the displacement are found in the upper 
and lower panels for the phonons in the diagonal and 
off-diagonal coupling baths, corresponding to the “polar¬ 
ized bath” and “unpolarized bath,” respectively. In the 
low frequency regimes, all the displacement coefficients 
fn,i and g n j converge to a value independent of n, i.e., 
fn,i = g-n .,1 —> c\i/2u)i in the polarized bath and 0 in the 
unpolarized bath, where c = —0.91 is a u >-independent 
constant. In the low frequency regime, however, f n j and 
g n j exhibit quite different behaviors, and the relations 
fi,i = = ~h,i and 9i,l = -92,1,92,1 = -93,1 are 



Figure 20: (Color on-line) The phase diagram of two-bath 
model is displayed in the as plane in the case of s = s and 
a = f). The phase boundary separating the critical phase 
from the localized phase is obtained from the parity index 
P(a, s) and the symmetry parameter of the rotational sym¬ 
metry 7 P h(a, s). 


found, indicating that the phonons in the unpolarized 
bath obey some kind of symmetry constraints. More¬ 
over, the quantum fluctuations of the two-bath model in 
the localized phase are negligible, since the amplitude of 
the high-frequency displacements A p « 0.6 in the unpo¬ 
larized bath is much smaller than that of the classical 
displacement |A;/2u;;| ~ 15 in the low-frequency limit. 

Fig. [22] shows the displacement coefficients f n g and 



Figure 21: (Color on-line) The displacement coefficients f n ,i 
and g n ,i of two-bath model in the localized phase are plotted 
as a function of the phonon frequency un at s = s = 0.5 
and a = (3 = 0.9. The upper and lower panels correspond to 
the polarized and unpolarized baths, respectively. The dashed 
line represents the classical displacements rescaled by a factor 
c = 0.91 for comparison. 
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Figure 22: (Color on-line) In the upper and lower panels, the 
displacement coefficients f n ,i and in the critical phase are 
plotted at s = s = 0.85, a = P = 0.02 for the polarized and 
unpolarized baths, respectively. The dashed line represents 
the classical displacements rescaled by a factor 0.87 for com¬ 
parison. 


g n ^i in the the critical phase as functions of u>i. Similar 
to the results in the localized phase shown in Fig. [21] 
one bath of the model is in the polarized state, and the 
other in the unpolarized state. However, the classical 
displacement f n j = g n j —> Xi/2uji « 0.2 at ug = 10 -6 
in the upper panel is compatible with the amplitude of 
the high-frequency displacements A p ps 0.15 in the lower 
panel. It means that the quantum fluctuations play an 
important role in the critical phase, unlike the case in the 
localized phase. 

Finally, the ground states of the two-bath model in the 
localized and critical phases are compared via the phonon 
population P p h(x, z) as shown in Fig. [531 According to 
the results in Figs. [21] and [221 the displacement coef¬ 
ficients f rlt i and g n , t i in the localized and critical phases 
are quite different in the low-frequency regime, especially 
at oji = 10 -6 . Hence, only the bath modes at the fre¬ 
quency uji = 10 -6 are considered, and the unit of the 
length 1 /yfuTi = 10 3 is set. The phonon state in the case 
of s = s = 0.85 and a = /3 = 0.02 (critical phase) is 
located nearby the origin O, but the one in the case of 
s = s = 0.5 and a = p = 0.9 (localized phase) is far 
away it. In both the cases, the distance d between the 
center of the phonon state and the origin is proportional 
to the displacement f n .i- g n ,i in the polarized bath. More¬ 
over, the central angle to the origin is calculated to be 
0 = 2arctan(r/d) ~ 0.04-7T for the case of the localized 
phase, where r = 0.4 is the radius of the phonon popula¬ 
tion P p h(x, z) and d = 20 is the distance. Interestingly, 
the central angle 0 is comparable with the width of the 
peaks A 9 = 0.0387T defined in Fig. [15] It further sup¬ 
ports that the parity index P = 2/S.O/ir < 1 reflects the 



x 

Figure 23: (Color on-line) The wave function of the ground 
states is displayed for the two cases of s = s = 0.85, a = 
P = 0.02 (nearby the origin O) and s = s = 0.5, a = /3 = 
0.9 (in a corner), corresponding to the critical and localized 
phases, respectively. The x- and 2-coordinate correspond to 
off-diagonal and diagonal coupling baths, respectively, and 
the colour represents the phonon population P p h(x,2) at the 
frequency uji = 10 -6 . For convenience, we set the unit of the 
length in the X-Z plane as 1/^/wj = 10 3 . 

localized nature of the ground state. 

B. The case with 

In the two bath model, we next investigate the case 
with (3 ^ a to further identify the critical and localized 
phases. According to Ref. [221, there exists a first-order 
quantum phase transition separating the doubly degen¬ 
erate “localized state” with |(< 7 Z )| > 0 and (a x ) = 0 from 
the doubly degenerate “delocalized state” with (a z ) = 0 
and |(<T a; )| > 0. The transition point /3 C = a is expected 
from the X-Z symmetry when the spectral exponents 
obey s = s. In the following variational calculations, 
we use N = 4 and M = 20, which have been found to be 
sufficient in obtaining reliable results of the localized-to- 
delocalized phase transition [22|] . 

In Fig. l24T a). the 2 component of the spin polariza¬ 
tion (cr 2 ) obtained by the variational method is plotted 
with respect to the ratio j3/a for various values of the 
diagonal coupling a in the case of s = s = 0.4. The tran¬ 
sition point is determined at /3 c /a = 1.0000(1), in perfect 
agreement with the expectation /3 C = a within numerical 
errors. Furthermore, a linear decay of ( a z ) is found for 
P < /3 C before showing an abrupt jump to zero at the 
transition point, thereby verifying the transition to be 
of the first order. For comparison, the numerical results 
obtained by the DMRG algorithm are also displayed in 
Fig. 135T b). Similar behavior of (a z ) is found, pointing to 
the validity of the variational method. 

Due to infinite degenerate ground states T(0)|4' s ), the 
spin polarization (a z ) at the transition point /3 C = a is 
variable as shown in Fig. 1131 In the regime with ft < a, 
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Figure 24: (Color on-line) The z component of the spin polar¬ 
ization (a z ) for various values of a is plotted as a function of 
the ratio /3/a in (a) and (b), corresponding to the variational 
and DMRG results, respectively. In both (a) and (b), the 
dash-dotted line indicates the transition point /3 C = a, and 
dashed lines represent linear fits. The value of the spectral 
exponent s = s = 0.4 is set. 
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however, the value of (a z ) is quite robust, since the ro¬ 
tational symmetry is broken by the anisotropic coupling. 
The generalized susceptibility y is thus introduced by 


X = 


dm 


T=0 


(34) 


where m is the spin polarization defined in Eq. ©, and 
r = |/3 — (3c\//3 c = |/3/a — 1| is the reduced coupling 
strength. As shown in Fig. [24j the £ component of the 
spin polarization (er 2 ) can be fitted by a linear behavior. 
Since the x component of the spin polarization is ( a x ) = 
0 in the localized state [22} , one obtains m = (<r 2 ) = 
or + b. The generalized susceptibility y = a can then be 
calculated with the linear fitting procedure for different 
values of a and s. An extended scaling form m = 1 — 
aexp(— br) could lead to a better fitting of the numerical 
data, and yield y = ab in the limit of r —>■ 0. Fig. l25l 
shows the generalized susceptibility y(a) at s = s = 0.4 
obtained by the variational method (solid triangles) and 
DMRG algorithm (open circles). Both of them decay 
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Figure 25: The generalized susceptibility y = dm/dr | 0 

is displayed as a function of the coupling strength a at 
s = s = 0.4. The circles and triangles correspond to the 
DMRG and variational results, respectively, and the dashed 
lines represent exponential fits. 


exponentially with a , though the exponents obtained via 
the fitting with y ~ exp(— ca d ) differ for the DMRG 
(d = 0.5) and variational method (d = 0.7). The two 
curves y(a) follow the same behavior, implying that the 
system is always in the localized phase at s = s = 0.4 for 
various values of a, and no phase transition occurs when 
s = 0.4 < s*. 

In Fig. ECTa). the generalized susceptibility y(a, s) is 
displayed as a function of the ratio s/(l + a) for various 
values of a. Unlike the trend shown in Fig. 1251 y(a, s) in¬ 
creases with the spectral exponent s until s/ (l+o) = 0.53 
marked by the dash-dotted line as the position of the 
peaks when a > 0.1. It points to the transition bound¬ 
ary separating the localized phase from the critical phase. 
Moreover, the generalized susceptibility y(a, s ) is also 
displayed in Fig. EBI b) for various values of s. The po¬ 
sition of the peak similarly marked by the dash-dotted 
line is found to be at s/(l + a) = 0.48 when s > 0.5. 
From these numerical results, one can obtain a relation¬ 
ship s/ (1 + cc) ~ 0.5, pointing to the transition boundary, 
i.e., a c = 2(s — 0.5). It further supports our contention 
that the critical value of the spectral exponent is s* ~ 0.5. 


V. CONCLUSION 

The ground states of the spin-boson model with di¬ 
agonal and off-diagonal coupling to two identical inde¬ 
pendent baths have been comprehensively studied in this 
paper by the variational approach, the DMRG algorithm 
and the exact diagonalization method. Adopting a gen¬ 
eralized trial wave function, i.e., multi-Di ansatz, as 
the variational wave function, the spin polarization m, 
ground state energy E g and wave function |4' 9 ) are cal- 
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Figure 26: (Color on-line) The generalized susceptibility 

X'(a, s) is displayed as a function of s/(l + a) in (a) for vari¬ 
ous values of a and in (b) for various values of s = 0.4. The 
dasli-dotted lines indicate the peak positions of the curves 
for a > 0.1 and s > 0.5, pointing to a transition boundary 
s/(l + a) ss 0.5 in the phase diagram. 


culated accurately by the variational method, in good 
agreement with those from the exact diagonalization in 
the case involving two oscillators and those from the 
DMRG algorithm in the case involving two baths de¬ 
scribing a continuous spectral density function. 

Three phases (localized, critical and free) are identi¬ 
fied, corresponding to the strong, intermediate and weak 
coupling regimes, respectively. Via the symmetry param¬ 
eters £ and 7 P h, the nature of the these phases is uncov¬ 
ered. The breakdown of the parity and rotational sym¬ 
metries is found to occur along the quantum phase transi¬ 
tion between the localized and critical phases. The phase 
boundary is determined by the parity index P(a, s) of the 
parity-symmetry regime with £ > 0, consistent with that 
obtained by 7 P h(a, s). Moreover, the critical value of the 
spectral exponent is estimated as s* = 0.49(1), well in 
agreement with the mean-field prediction 1/2 [22[. The 
behavior of the spin polarization m(a, s) and the general¬ 
ized susceptibility x( a > s ) is also investigated for various 
values of the coupling strengthes a and spectral expo¬ 
nents s. Both the results point to s* ~ 0.5, further sup¬ 
porting the accuracy of our results in pinning the transi¬ 
tion point. 
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